ggplot geospatial example

Geospartial Data Visualization with ggplot

Author

sungil_park

Published

March 15, 2023

Modified

March 20, 2023

library packages

rm(list = ls())
library(ggplot2)
library(dplyr)


다음의 패키지를 부착합니다: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(colorspace)
library(stringr)
library(geojsonsf)
library(sf)

Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE

Map Korea

load `.json`

kor_sigu <- geojson_sf("https://raw.githubusercontent.com/Sungileo/trainsets/main/kor/KOR_SIGU.json")

Warning in readLines(con):
'https://raw.githubusercontent.com/Sungileo/trainsets/main/kor/KOR_SIGU.json'에서
불완전한 마지막 행이 발견되었습니다

kor_202202 <- read.csv("https://raw.githubusercontent.com/Sungileo/trainsets/main/202202_%EC%A3%BC%EB%AF%BC%EB%93%B1%EB%A1%9D%EC%9D%B8%EA%B5%AC%EB%B0%8F%EC%84%B8%EB%8C%80%ED%98%84%ED%99%A9.csv")

kor_202202 %>% head()

           행정구역 행정구역_코드 총인구수  세대수 세대당_인구 남자_인구수
1       서울특별시     1100000000  9508451 4442586        2.14     4615823
2 서울특별시 종로구    1111000000   144575   73763        1.96       70092
3   서울특별시 중구    1114000000   122167   63644        1.92       59446
4 서울특별시 용산구    1117000000   222413  111134        2.00      106881
5 서울특별시 성동구    1120000000   285137  134286        2.12      138866
6 서울특별시 광진구    1121500000   340494  168975        2.02      164226
  여자_인구수 남여_비율
1     4892628      0.94
2       74483      0.94
3       62721      0.95
4      115532      0.93
5      146271      0.95
6      176268      0.93

kor_202202 %>% sapply(class)

     행정구역 행정구역_코드      총인구수        세대수   세대당_인구 
  "character"     "numeric"     "numeric"     "numeric"     "numeric" 
  남자_인구수   여자_인구수     남여_비율 
    "numeric"     "numeric"     "numeric"

kor_202202$행정구역_코드 <- kor_202202$행정구역_코드 %>% format()

Merge data by 행정구역_코드

use_map <- kor_sigu
use_map$행정구역_코드 <- paste(use_map$SIG_CD,"00000",sep = "")
use_map <- use_map %>% merge(kor_202202,by = "행정구역_코드", all.x=T)

Plot

use_map %>% ggplot(aes(fill=총인구수))+
  geom_sf(color = "grey90")+
  coord_sf(datum = NA)+
  scale_fill_distiller(
    name = "인구수",
    palette = "Blues", type = "seq", na.value = "grey60",
    direction = 1,
    breaks = seq(0,10,2) * 1e+5,
    labels = format(seq(0,10,2) * 1e+5, big.mark = ",",scientific = FALSE))+
  theme_minimal()+
  theme(
    legend.title.align = 0.5,
    legend.text.align = 1.0,
    legend.position = c(0.85,0.2)
  )

Plot Daejeon

daejeon_map <-  use_map %>% filter(행정구역 %>% substr(1,5) == "대전광역시")

daejeon_map %>% ggplot(aes(fill=총인구수))+
  geom_sf(color = "grey90")+
  coord_sf(datum = NA)+
  scale_fill_distiller(
    name = "인구수",
    palette = "Blues", type = "seq", na.value = "grey60",
    direction = 1,
    breaks = seq(0,10,2) * 1e+5,
    labels = format(seq(0,10,2) * 1e+5, big.mark = ",",scientific = FALSE))+
  theme_minimal()+
  theme(
    legend.title.align = 0.5,
    legend.text.align = 1.0,
    legend.position = c(0.95,0.2)
  )

Plot Gender_ratio

use_map %>% ggplot(aes(fill = 남여_비율))+
  geom_sf()+
  scale_fill_continuous_diverging(
    name = "남자/여자",
    palette = "BLue-Red",
    mid=1,
  limits = 1 + c(-1,+1)*0.35,
  rev = T)+
  theme_minimal()+
  theme(
    legend.title.align = 0.5,
    legend.text.align = 1.0,
    legend.position = c(0.85,0.2)
  )

2023년 3월 시군구 인구수

시군구 지도 데이터, 행정구역 코드 10자리로 만들기

use_map <- kor_sigu
use_map$행정구역_코드 <- paste(use_map$SIG_CD,"00000",sep = "")

데이터 로드

2023년 3월 주민등록 인구통계 데이터, 행정안전부

data_pop <- read.csv("https://raw.githubusercontent.com/Sungileo/trainsets/main/202303_202303_%EC%A3%BC%EB%AF%BC%EB%93%B1%EB%A1%9D%EC%9D%B8%EA%B5%AC%EB%B0%8F%EC%84%B8%EB%8C%80%ED%98%84%ED%99%A9_%EC%9B%94%EA%B0%84.csv",encoding = "utf-8")
data_pop %>% head()

                        행정구역 X2023년03월_총인구수 X2023년03월_세대수
1       서울특별시  (1100000000)            9,426,404          4,463,385
2 서울특별시 종로구 (1111000000)              141,060             72,679
3   서울특별시 중구 (1114000000)              120,963             63,862
4 서울특별시 용산구 (1117000000)              217,756            109,735
5 서울특별시 성동구 (1120000000)              280,240            133,513
6 서울특별시 광진구 (1121500000)              336,801            169,787
  X2023년03월_세대당.인구 X2023년03월_남자.인구수 X2023년03월_여자.인구수
1                    2.11               4,566,299               4,860,105
2                    1.94                  68,170                  72,890
3                    1.89                  58,699                  62,264
4                    1.98                 104,640                 113,116
5                    2.10                 136,233                 144,007
6                    1.98                 162,209                 174,592
  X2023년03월_남여.비율
1                  0.94
2                  0.94
3                  0.94
4                  0.93
5                  0.95
6                  0.93

전처리

인구수0명 출장소 제외
행정구역 코드 10자리 추출
인구수 숫자형으로 변환
시 단위 제외, 정렬

data_202303 <- data_pop %>% 
  filter(X2023년03월_총인구수>0) %>%  
  select(행정구역,X2023년03월_총인구수) %>% 
  mutate(행정구역_코드 = str_sub(행정구역,-11,-2),
         X2023년03월_총인구수 = gsub(",","",X2023년03월_총인구수) %>% as.numeric()) %>% 
  filter(substr(행정구역_코드,3,4)!="00") %>% 
  arrange(desc(X2023년03월_총인구수))

data_202303 %>% head()

                      행정구역 X2023년03월_총인구수 행정구역_코드
1   경기도 수원시 (4111000000)              1192225    4111000000
2   경기도 고양시 (4128000000)              1077909    4128000000
3   경기도 용인시 (4146000000)              1073513    4146000000
4 경상남도 창원시 (4812000000)              1017273    4812000000
5   경기도 성남시 (4113000000)               923416    4113000000
6   경기도 화성시 (4159000000)               922231    4159000000

지도 데이터와 병합

use_map <- use_map %>% merge(data_202303,by = "행정구역_코드", all.x=T)

Plot

use_map %>% ggplot(aes(fill=X2023년03월_총인구수))+
  geom_sf(color = "grey90")+
  coord_sf(datum = NA)+
  scale_fill_distiller(
    name = "2023년 3월 인구수",
    palette = "Blues", type = "seq", na.value = "grey60",
    direction = 1,
    breaks = seq(0,10,2) * 1e+5,
    labels = format(seq(0,10,2) * 1e+5, big.mark = ",",scientific = FALSE))+
  theme_minimal()+
  theme(
    legend.title.align = 0.5,
    legend.text.align = 1.0,
    legend.position = c(0.85,0.2))+
  labs(title = "2023년 3월")

인구수 증가율

file_2023 <- read.csv("https://raw.githubusercontent.com/Sungileo/trainsets/main/202303_202303_%EC%A3%BC%EB%AF%BC%EB%93%B1%EB%A1%9D%EC%9D%B8%EA%B5%AC%EB%B0%8F%EC%84%B8%EB%8C%80%ED%98%84%ED%99%A9_%EC%9B%94%EA%B0%84.csv",encoding = "utf-8")

file_2013 <- read.csv("https://raw.githubusercontent.com/Sungileo/trainsets/main/201303_201303_%EC%A3%BC%EB%AF%BC%EB%93%B1%EB%A1%9D%EC%9D%B8%EA%B5%AC%EB%B0%8F%EC%84%B8%EB%8C%80%ED%98%84%ED%99%A9_%EC%9B%94%EA%B0%84.csv",encoding = "UTF-8")

전처리

인구수 0명 이상 필터
행정구역(지역코드), 총인구수 열만 선택
행정구역(지역코드)에서 지역코드와 행정구역 분리
시 단위 제외, 인구수 기준 정렬

data_2023 <- file_2023 %>%
  filter(X2023년03월_총인구수>0) %>%  
  select(행정구역,X2023년03월_총인구수) %>% 
  mutate(행정구역_코드 = str_sub(행정구역,-11,-2),
         X2023년03월_총인구수 = gsub(",","",X2023년03월_총인구수) %>% as.numeric(),
         행정구역 =  sapply(행정구역, function(x) strsplit(x, "(", fixed=T)[[1]][1]),
         행정구역 = sapply(행정구역, function(x) gsub("( *)$", "", x) %>% paste())) %>% 
  filter(substr(행정구역_코드,3,4)!="00") %>% 
  arrange(desc(X2023년03월_총인구수))

data_2013 <- file_2013 %>% 
  filter(X2013년03월_총인구수>0) %>%  
  select(행정구역,X2013년03월_총인구수) %>% 
  mutate(행정구역_코드 = str_sub(행정구역,-11,-2),
         X2013년03월_총인구수 = gsub(",","",X2013년03월_총인구수) %>% as.numeric(),
         행정구역 =  sapply(행정구역, function(x) strsplit(x, "(", fixed=T)[[1]][1]),
         행정구역 = sapply(행정구역, function(x) gsub("( *)$", "", x) %>% paste())) %>% 
  filter(substr(행정구역_코드,3,4)!="00") %>% 
  arrange(desc(X2013년03월_총인구수))

병합

지역코드 기준 병합
인구성장률 열 추가
중복 열 제거, 인구성장률 기준 정렬
서울, 대전, 대구, 부산지역만 필터
시도 추출, factor 변환

kor_census <- data_2013 %>% 
  merge(data_2023,by = "행정구역_코드", all.x=T) %>%  
  mutate(성장률 = (X2023년03월_총인구수 - X2013년03월_총인구수) / X2013년03월_총인구수) %>% 
  select(행정구역.x,X2013년03월_총인구수,X2023년03월_총인구수, 성장률, 행정구역_코드) %>% 
  filter(substr(행정구역.x,1,2) %in% c("서울","대전","대구","부산")) %>%
  arrange(desc(성장률))

names(kor_census) <- c("행정구역", "X2013인구수","X2023인구수","성장률","행정구역_코드")

kor_census$시도 = sapply(kor_census$행정구역,
                           function(x) strsplit(x, " ")[[1]][1])
kor_census$시도 = factor(kor_census$시도,
                           levels = c("서울특별시","대전광역시","대구광역시","부산광역시"))

Plot

region_colors <- c("#E69F00","#56B4E9","#009E73","#F0E442")


ggplot(kor_census,aes(x = reorder(행정구역,성장률),y= 성장률, fill = 시도))+
  geom_col()+
  scale_y_continuous(name = "인구성장률",
                     expand = c(0,0),
                     labels = scales::percent_format(scale = 100))+
  scale_fill_manual(values = region_colors)+
  coord_flip()+
  theme_light()+
  theme(panel.border = element_blank(),
        panel.grid.major.y = element_blank())+
  theme(axis.title.y = element_blank(),
        axis.line.y = element_blank(),
        axis.ticks.length = unit(0,"pt"),
        axis.text.y = element_text(size = 8),legend.position = c(.78,.28),legend.background = element_rect(fill = "#FFFFFFB0"))

Map plot

kor_map <- kor_sigu
kor_map$행정구역_코드 <- paste(kor_map$SIG_CD,"00000",sep="")

kor_census_2 <- data_2013 %>% 
  merge(data_2023,by = "행정구역_코드", all.x=T) %>%  
  mutate(성장률 = (X2023년03월_총인구수 - X2013년03월_총인구수) / X2013년03월_총인구수) %>% 
  select(행정구역.x,X2013년03월_총인구수,X2023년03월_총인구수, 성장률, 행정구역_코드) %>% 
  arrange(desc(성장률))

names(kor_census_2) <- c("행정구역", "X2013인구수","X2023인구수","성장률","행정구역_코드")


kor_map <- kor_map %>% left_join(kor_census_2, by="행정구역_코드")

kor_map %>% ggplot(aes(fill=성장률))+
  geom_sf()+
  scale_fill_continuous_diverging(
    name = "인구성장률",
    palette = "BLue-Red",
    limits = c(-0.4,2.4))+
  theme_minimal()+
  theme(legend.title.align = 0.5,
        legend.text.align = 1.0,
        legend.position = c(0.85,0.2))

인구증감

data <- read.csv("https://raw.githubusercontent.com/Sungileo/trainsets/main/kor_census_2013_2023.csv",encoding = "utf-8")

data <- data %>% mutate(인구증감 = 총인구수_2023-총인구수_2013)

use_map <- kor_sigu
use_map$행정구역_코드 <- paste(use_map$SIG_CD,"00000",sep = "") %>% as.numeric()
use_map <- use_map %>% merge(data,by = "행정구역_코드")

use_map %>% ggplot(aes(fill = 인구증감))+
  geom_sf()+
  coord_sf(datum = NA)+
  scale_fill_continuous_diverging(
    name = "인구증감",
    palette = "BLue-Red",
    na.value = "grey40",
    mid=0,
    rev = T,
    limits = c(-4,4)*100000,
    labels = format(seq(-4,4,2) * 1e+5, big.mark = ",",scientific = FALSE))+
  theme_minimal()+
  theme(legend.position = c(0.85,0.2))

Mapdeck

library(dplyr)
library(mapdeck)
set_token("pk.eyJ1Ijoic3VuZ2lsZW8iLCJhIjoiY2xoYTRwbXEzMGR6eTNkbXBoZnluNXdyYSJ9.Id1fKIbhtvA9Mrnyo_1JQA")

crash_data = read.csv("https://git.io/geocompr-mapdeck")
crash_data = na.omit(crash_data)

ms = mapdeck_style("dark")

# mapdeck(style = "mapbox://styles/sungileo/clhtwpzyw00rl01rhfoiafw6q")

grid_map

mapdeck()함수의 주요 Arguments
- style : the style of the map (see mapdeck_style) : one of streets, outdoors, light, dark, satellite, satellite-streets 생성할 맵의 스타일
- pitch : the pitch angle of the map 맵의 기울기를 설정합니다. 0은 수평을 의미하고, 60은 수직
- zoom : zoom level of the map 초기 맵의 줌 레벨을 설정
- location : unnamed vector of lon and lat coordinates (in that order) 초기 맵의 중심 위치를 지정
add_grid()함수의 주요 Argurments
- data: 격자를 생성하는 데 사용할 데이터 프레임
- lon : column containing longitude values
- lat : column containing latitude values
- cell_size : size of each cell in meters. Default 1000 격자 셀의 크기
- elevation : the height the polygon extrudes from the map. Only available if neither stroke_colour or stroke_width are supplied. Default 0 격자의 높이
- layer_id : single value specifying an id for the layer 격자 레이어의 고유 식별자를 지정하는 데 사용. 이 식별자는 나중에 해당 레이어를 참조하거나 다른 함수를 통해 조작할 때 사용
- colour_range : vector of 6 hex colours

mapdeck(style = ms, 
        pitch = 45, 
        location = c(0, 52), 
        zoom = 4) %>%
  add_grid(data = crash_data, 
           lat = "lat", 
           lon = "lng", 
           cell_size = 1000,
           elevation_scale = 50, 
           layer_id = "grid_layer",
           colour_range = viridisLite::plasma(6))

Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite

arc map

add_arc()함수의 주요 Argurments

add_arc function - RDocumentation
- layer_id : single value specifying an id for the layer. Use this value to distinguish between shape layers of the same type. Layers with the same id are likely to conflict and not plot correctly
- origin : vector of longitude and latitude columns, and optionally an elevation column, or an sfc column
- destination : vector of longitude and latitude columns, and optionally an elevatino column, or an sfc column
- stroke_from : column of data or hex colour to use as the staring stroke colour. IIf using a hex colour, use either a single value, or a column of hex colours on data
- stroke_to : column of data or hex colour to use as the ending stroke colour. If using a hex colour, use either a single value, or a column of hex colours on data
- stroke_width : width of the stroke in pixels

Visualize flight data in USA

url <- 'https://raw.githubusercontent.com/plotly/datasets/master/2011_february_aa_flight_paths.csv'
flights <- read.csv(url)
flights$id <- seq_len(nrow(flights))
flights$stroke <- sample(1:3, size = nrow(flights), replace = T)

mapdeck(style = ms, 
        pitch = 45 ) %>%
  add_arc(data = flights, 
          layer_id = "arc_layer", 
          origin = c("start_lon", "start_lat"), 
          destination = c("end_lon", "end_lat"), 
          stroke_from = "stroke", 
          stroke_to = "stroke", 
          stroke_width = "stroke" 
  )

mapdeck(style = 'mapbox://styles/mapbox/dark-v9', 
        pitch = 45 ) %>%
  add_animated_arc(data = flights, 
                   layer_id = "arc_layer", 
                   origin = c("start_lon", "start_lat"), 
                   destination = c("end_lon", "end_lat"), 
                   stroke_from = "stroke", 
                   stroke_to = "stroke", 
                   stroke_width = "stroke")

animated_arc is an experimental layer and the function may change without warning

heatmap

mapdeck(style = ms, 
        pitch = 45, 
        location = c(0, 52), 
        zoom = 4)%>%
  add_heatmap(data = crash_data[1:30000,],
              lat = "lat",
              lon = "lng", 
              colour_range = colourvalues::colour_values(1:6, palette = "inferno"))

roads %>% View()

mapdeck(style = mapdeck_style("dark"), 
        zoom = 10) %>%
  add_path(data = roads,
            stroke_colour = "FQID",
           layer_id = "path_layer")

hexagon map

mapdeck(style = mapdeck_style("dark"), 
        pitch = 45, 
        location = c(0, 52), 
        zoom = 4) %>%
  add_hexagon(data = crash_data, 
           lat = "lat", 
           lon = "lng", 
           elevation_scale = 100, 
           layer_id = "hex_layer",
           colour_range = viridisLite::plasma(6))

leaflet

library(leaflet)
library(spData)

The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
     (status 2 uses the sf package in place of rgdal)

cycle_hire %>% head()

Simple feature collection with 6 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -0.1975742 ymin: 51.49313 xmax: -0.08460569 ymax: 51.53006
Geodetic CRS:  WGS 84
  id               name             area nbikes nempty
1  1       River Street      Clerkenwell      4     14
2  2 Phillimore Gardens       Kensington      2     34
3  3 Christopher Street Liverpool Street      0     32
4  4  St. Chad's Street     King's Cross      4     19
5  5     Sedding Street    Sloane Square     15     12
6  6 Broadcasting House       Marylebone      0     18
                      geometry
1  POINT (-0.1099705 51.52916)
2  POINT (-0.1975742 51.49961)
3 POINT (-0.08460569 51.52128)
4  POINT (-0.1209737 51.53006)
5   POINT (-0.156876 51.49313)
6  POINT (-0.1442289 51.51812)

lnd %>% head()

Simple feature collection with 6 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -0.4615023 ymin: 51.28676 xmax: 0.3340155 ymax: 51.63173
Geodetic CRS:  WGS 84
                  NAME  GSS_CODE  HECTARES NONLD_AREA ONS_INNER SUB_2009
1 Kingston upon Thames E09000021  3726.117      0.000         F     <NA>
2              Croydon E09000008  8649.441      0.000         F     <NA>
3              Bromley E09000006 15013.487      0.000         F     <NA>
4             Hounslow E09000018  5658.541     60.755         F     <NA>
5               Ealing E09000009  5554.428      0.000         F     <NA>
6             Havering E09000016 11445.735    210.763         F     <NA>
  SUB_2006                       geometry
1     <NA> MULTIPOLYGON (((-0.3306791 ...
2     <NA> MULTIPOLYGON (((-0.06402124...
3     <NA> MULTIPOLYGON (((0.01213094 ...
4     <NA> MULTIPOLYGON (((-0.2445624 ...
5     <NA> MULTIPOLYGON (((-0.4118327 ...
6     <NA> MULTIPOLYGON (((0.1586928 5...

pal = colorNumeric("RdYlBu", domain = cycle_hire$nbikes)

leaflet(data = cycle_hire) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircles(col = ~pal(nbikes), opacity = 0.9) %>% 
  addPolygons(data = lnd, fill = FALSE) %>% 
  addLegend(pal = pal, values = ~nbikes) %>% 
  setView(lng = -0.1, 51.5, zoom = 12) %>% 
  addMiniMap()

# creat a basic map
leaflet() %>%
  addTiles() %>% # add default OpenStreetMap map tiles
  setView( lng=127.063, lat=37.513, zoom = 6) # korea, zoom 6

지도 스타일 추가하기
- addProviderTiles(“Tile Name Here”) 를 이용하여 외부 지도 타일을 추가(아래에서 원하는 스타일 map 이름 선택)

# map style: NASA
leaflet() %>%
  addTiles() %>%
  setView( lng=127.063, lat=37.513, zoom = 6) %>%
  addProviderTiles("NASAGIBS.ViirsEarthAtNight2012")

# map style: Esri.WorldImagery
leaflet() %>%
  addTiles() %>%
  setView( lng=127.063, lat=37.513, zoom = 16) %>%
  addProviderTiles("Esri.WorldImagery")

지도 좌표에 표시 풍선 추가하기
- addMarkers() 메소드를 사용하면 위도(latitude), 경도(longitude) 좌표 위치에 풍선 모양의 표식과 커서를 클릭했을 때 팝업으로 나타나는 설명을 추가할 수 있음

# adding Popup
popup = c("한남대학교 빅데이터응용학과")

leaflet() %>%
  addTiles() %>%
  addMarkers(lng = c(127.4219), # longitude
             lat = c(36.3548), # latitude
             popup = popup)

ggmap

library(ggmap)

ℹ Google's Terms of Service: <https://mapsplatform.google.com>
ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.

getmap <- get_googlemap("seoul")

ℹ <https://maps.googleapis.com/maps/api/staticmap?center=seoul&zoom=10&size=640x640&scale=2&maptype=terrain&key=xxx>

ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=seoul&key=xxx>

ggmap(getmap)

str(crime)

'data.frame':   86314 obs. of  17 variables:
 $ time    : POSIXct, format: "2010-01-01 15:00:00" "2010-01-01 15:00:00" ...
 $ date    : chr  "1/1/2010" "1/1/2010" "1/1/2010" "1/1/2010" ...
 $ hour    : int  0 0 0 0 0 0 0 0 0 0 ...
 $ premise : chr  "18A" "13R" "20R" "20R" ...
 $ offense : Factor w/ 7 levels "aggravated assault",..: 4 6 1 1 1 3 3 3 3 3 ...
 $ beat    : chr  "15E30" "13D10" "16E20" "2A30" ...
 $ block   : chr  "9600-9699" "4700-4799" "5000-5099" "1000-1099" ...
 $ street  : chr  "marlive" "telephone" "wickview" "ashland" ...
 $ type    : chr  "ln" "rd" "ln" "st" ...
 $ suffix  : chr  "-" "-" "-" "-" ...
 $ number  : int  1 1 1 1 1 1 1 1 1 1 ...
 $ month   : Ord.factor w/ 8 levels "january"<"february"<..: 1 1 1 1 1 1 1 1 1 1 ...
 $ day     : Ord.factor w/ 7 levels "monday"<"tuesday"<..: 5 5 5 5 5 5 5 5 5 5 ...
 $ location: chr  "apartment parking lot" "road / street / sidewalk" "residence / house" "residence / house" ...
 $ address : chr  "9650 marlive ln" "4750 telephone rd" "5050 wickview ln" "1050 ashland st" ...
 $ lon     : num  -95.4 -95.3 -95.5 -95.4 -95.4 ...
 $ lat     : num  29.7 29.7 29.6 29.8 29.7 ...

Houstonmap <- get_map("Houston")

ℹ <https://maps.googleapis.com/maps/api/staticmap?center=Houston&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx>

ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Houston&key=xxx>

ggmap(Houstonmap)

ggmap(Houstonmap) + geom_point(data=crime, aes(x=lon,y=lat), size=0.1, alpha=0.1) #점의크기, 점의 투명도 조절

Warning: Removed 237 rows containing missing values (`geom_point()`).

#지도 확대 & 특정 지역 데이터만 추출하기
Houstonmap <- get_map("Houston", zoom=14)

ℹ <https://maps.googleapis.com/maps/api/staticmap?center=Houston&zoom=14&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx>

ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Houston&key=xxx>

crime1 <- crime[(crime$lon < -95.344 & crime$lon > -95.395) & 
                  (crime$lat < 29.783 & crime$lat > 29.738),]
ggmap(Houstonmap) + geom_point(data=crime1, aes(x=lon,y=lat), alpha=0.3)

Warning: Removed 5 rows containing missing values (`geom_point()`).

ggmap(Houstonmap) + geom_point(data=crime1, aes(x=lon, y=lat, colour = offense))

Warning: Removed 5 rows containing missing values (`geom_point()`).

crime2 <- crime1[!duplicated(crime1[,c("lon","lat")]),] #위,경도에 대해 중복되지 않게 하나의 관측치만 선택
crime2$offense <- as.character(crime2$offense) #범죄 종류 문자형으로 변경

crime2$offense[crime2$offense=="murder" | crime2$offense=="rape"] <- "4"
crime2$offense[crime2$offense=="robbery" | crime2$offense=="aggravated assault"] <- "3"
crime2$offense[crime2$offense=="burglary" | crime2$offense=="auto theft"] <- "2"
crime2$offense[crime2$offense=="theft"] <- "1"

crime2$offense <-  as.numeric(crime2$offense) #문자형을 숫자로 변경

ggmap(Houstonmap) + geom_point(data=crime2, aes(x=lon, y=lat, size = offense),alpha=0.2)

Warning: Removed 1 rows containing missing values (`geom_point()`).

#범죄 위험도에 따라 점의 크기 및 색깔로 구별
ggmap(Houstonmap) + geom_point(data=crime2, aes(x=lon, y=lat, size = offense, colour=offense),
                               alpha=0.5) +
  scale_colour_gradient(low="white", high="red")

Warning: Removed 1 rows containing missing values (`geom_point()`).

crime3 <- crime2[crime2$date=="1/1/2010",]
crime4 <- crime3[!duplicated(crime3[,c("hour")]),]
ggmap(Houstonmap) + geom_point(data=crime3, aes(x=lon, y=lat)) +
  geom_text(data=crime4, aes(label=street), vjust=1.2) +
  geom_path(data=crime4, aes(x=lon, y=lat), color="red")

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing missing values (`geom_text()`).

Warning: Removed 1 row containing missing values (`geom_path()`).

library(sf)
library(dplyr)
library(spDataLarge)
library(stplanr)      # geographic transport data package


다음의 패키지를 부착합니다: 'stplanr'

The following object is masked from 'package:ggmap':

    route

library(tmap)         # visualization package (see Chapter 8)

names(bristol_zones)

[1] "geo_code" "name"     "geometry"

bristol_zones

Simple feature collection with 102 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -2.845847 ymin: 51.28248 xmax: -2.252388 ymax: 51.73982
Geodetic CRS:  WGS 84
First 10 features:
      geo_code                             name                       geometry
2905 E02002985 Bath and North East Somerset 001 MULTIPOLYGON (((-2.510462 5...
2907 E02002987 Bath and North East Somerset 003 MULTIPOLYGON (((-2.476122 5...
2925 E02003005 Bath and North East Somerset 021 MULTIPOLYGON (((-2.55073 51...
2932 E02003012                      Bristol 001 MULTIPOLYGON (((-2.595763 5...
2933 E02003013                      Bristol 002 MULTIPOLYGON (((-2.593783 5...
2934 E02003014                      Bristol 003 MULTIPOLYGON (((-2.639581 5...
2935 E02003015                      Bristol 004 MULTIPOLYGON (((-2.584973 5...
2936 E02003016                      Bristol 005 MULTIPOLYGON (((-2.565948 5...
2937 E02003017                      Bristol 006 MULTIPOLYGON (((-2.616485 5...
2938 E02003018                      Bristol 007 MULTIPOLYGON (((-2.637681 5...

nrow(bristol_od)

[1] 2910

bristol_od

# A tibble: 2,910 × 7
   o         d           all bicycle  foot car_driver train
   <chr>     <chr>     <dbl>   <dbl> <dbl>      <dbl> <dbl>
 1 E02002985 E02002985   209       5   127         59     0
 2 E02002985 E02002987   121       7    35         62     0
 3 E02002985 E02003036    32       2     1         10     1
 4 E02002985 E02003043   141       1     2         56    17
 5 E02002985 E02003049    56       2     4         36     0
 6 E02002985 E02003054    42       4     0         21     0
 7 E02002985 E02003100    22       0     0         19     3
 8 E02002985 E02003106    48       3     1         33     8
 9 E02002985 E02003108    31       0     0         29     1
10 E02002985 E02003121    42       1     2         34     0
# ℹ 2,900 more rows

nrow(bristol_zones)

[1] 102

#bristol_zones 지역별로 출발지 'o'로 그룹화 하여 이동량 계산
zones_attr = bristol_od %>% 
  group_by(o) %>% 
  summarize_if(is.numeric, sum) %>% 
  dplyr::rename(geo_code = o)

summary(zones_attr$geo_code %in% bristol_zones$geo_code)

   Mode    TRUE 
logical     102

#bristol_zones과 zones_attr를  geo_code를 기준으로 합치기

zones_joined = left_join(bristol_zones, zones_attr, by = "geo_code")
sum(zones_joined$all)

[1] 238805

names(zones_joined)

[1] "geo_code"   "name"       "all"        "bicycle"    "foot"      
[6] "car_driver" "train"      "geometry"

names(zones_joined)[3] <- c("all_orig")
names(zones_joined)

[1] "geo_code"   "name"       "all_orig"   "bicycle"    "foot"      
[6] "car_driver" "train"      "geometry"

# 목적지가 geo_code인 데이터 group화하여 이동량 계산, 그리고 zones_joined에 추가하기
zones_od = bristol_od %>% 
  group_by(d) %>% 
  summarize_if(is.numeric, sum) %>% 
  dplyr::select(geo_code = d, all_dest = all) %>% 
  inner_join(zones_joined, ., by = "geo_code")

qtm(zones_od, c("all_orig", "all_dest")) +
  tm_layout(panel.labels = c("Origin", "Destination"))

od_top5 = bristol_od %>% 
  arrange(desc(all)) %>% 
  top_n(5, wt = all)

#활동성(전체 이동수단인원에서 자전거or 도보로 이동한 인원)
bristol_od$Active = (bristol_od$bicycle + bristol_od$foot) /
  bristol_od$all * 100

od_intra = filter(bristol_od, o == d) #o와 d가 같을때
od_inter = filter(bristol_od, o != d) #o와 d가 다를때

desire_lines = od2line(od_inter, zones_od)

Creating centroids representing desire line start and end points.

#> Creating centroids representing desire line start and end points.
qtm(desire_lines, lines.lwd = "all")

Legend labels were too wide. Therefore, legend.text.size has been set to 0.61. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

tmap_mode("plot")

tmap mode set to plotting

desire_lines_top5 = od2line(od_top5, zones_od)

Creating centroids representing desire line start and end points.

# tmaptools::palette_explorer()

tm_shape(desire_lines) +
  tm_lines(palette = "plasma", 
           breaks = c(0, 5, 10, 20, 40, 100),
           lwd = "all",
           scale = 9,
           title.lwd = "Number of trips",
           alpha = 0.6,
           col = "Active",
           title = "Active travel (%)") +
  tm_shape(desire_lines_top5) +
  tm_lines(lwd = 5, col = "black", alpha = 0.7) +
  tm_scale_bar()

Legend labels were too wide. Therefore, legend.text.size has been set to 0.61. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

desire_lines$distance = as.numeric(st_length(desire_lines))

desire_carshort = dplyr::filter(desire_lines, car_driver > 300 & distance < 5000)

route_carshort = route(l = desire_carshort, route_fun = route_osrm)

Most common output is sf

desire_carshort$geom_car = st_geometry(route_carshort)

plot(st_geometry(desire_carshort))
plot(desire_carshort$geom_car, col = "red", add = TRUE)
plot(st_geometry(st_centroid(zones_od)), add = TRUE)

Warning: st_centroid assumes attributes are constant over geometries

getmap <- get_googlemap("bristol", zoom = 11)

ℹ <https://maps.googleapis.com/maps/api/staticmap?center=bristol&zoom=11&size=640x640&scale=2&maptype=terrain&key=xxx>

ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=bristol&key=xxx>

bristol_map <- ggmap(getmap)
#센터 조정
getmap <- get_googlemap(center = c(-2.56, 51.53), zoom = 12)

ℹ <https://maps.googleapis.com/maps/api/staticmap?center=51.53,-2.56&zoom=12&size=640x640&scale=2&maptype=terrain&key=xxx>

bristol_map <- ggmap(getmap)

bristol_map + 
  geom_sf(data = desire_carshort, inherit.aes = F) + 
  geom_sf(data = desire_carshort$geom_car, inherit.aes = F, col="red") +
  geom_sf(data = st_geometry(st_centroid(zones_od)), inherit.aes = F)

Coordinate system already present. Adding new coordinate system, which will
replace the existing one.

Warning: st_centroid assumes attributes are constant over geometries

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

tmap_mode("view")

tmap mode set to interactive viewing

tm_basemap("OpenStreetMap")+
  tm_shape(desire_carshort$geometry)+
  tm_lines()+
  tm_shape(desire_carshort$geom_car)+
  tm_lines(col = "red",lwd = 0.5)+
  tm_shape(st_geometry(st_centroid(zones_od)))+
  tm_dots()

Warning: st_centroid assumes attributes are constant over geometries

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

Warning in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy, :
GDAL Error 1: PROJ: proj_as_wkt: DatumEnsemble can only be exported to
WKT2:2019

library packages

Map Korea

load .json

Merge data by 행정구역_코드

Plot

Plot Daejeon

Plot Gender_ratio

2023년 3월 시군구 인구수

데이터 로드

전처리

지도 데이터와 병합

Plot

인구수 증가율

전처리

병합

Plot

Map plot

인구증감

Mapdeck

grid_map

arc map

heatmap

hexagon map

leaflet

ggmap

load `.json`